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A genera iterative method proposed some years ago for the description of relativistic collapse, 
is presented here in comoving coordinates. For doing that we redefine the basic concepts required 
for the implementation of the method for comoving coordinates. In particular the definition of the 
post-quasistatic approximation in comoving coordinates is given. We write the field equations, the 
boundary conditions and a set of ordinary differential equations (the surface equations) which play 
a fundamental role in the algorithm. As an illustration of the method, we show how to build up 
a model inspired in the well known Schwarzschild interior solution. Both, the adiabatic and non 
adiabatic, cases are considered. 
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I. INTRODUCTION 

One of the most outstanding problems in the relativis- 
tic astrophysics and gravitation theory today is to pro- 
vide an accurate description of the gravitational collapse 
of a supermassive star. The final fate of such process [T]- 
[TO] (naked singularities, black holes, anything else ?), 
the mechanism behind a type II supernova event [llj- 
[TT) or the structure and evolution of the compact object 
resulting from such a process [THJ-dU], stand among the 
most interesting questions associated to that problem. 

There are essentially two approaches to describe the 
gravitational collapse in the context of general relativity. 
On the one hand, one may resort to numerical meth- 
ods which allow for considering more realistic 
equations of state. However, the obtained results, in gen- 
eral, are restricted and highly model dependents. Also, 
specific difficulties, associated to numerical solutions of 
partial differential equations in presence of shocks may 
complicate further the problem. It would be desirable 
in some cases less complicated numerical solvers that 
community could handle and adapt easily. On the other 
hand, one can use analytical solutions to Einstein equa- 
tions, which are more suitable for a general discussion, 
and may be very useful in the study of the structure and 
evolution of self-gravitating systems, since they may be 
relatively simple to analyze but still may contain some 
of the essential features of a realistic situation. However, 
often they are found, either for too simplistic equations 
of state and/or under additional heuristic assumptions 
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whose justification is usually uncertain. Occasionally, 
analytical approaches [2S]-[3n] challenge practical ones, 
say numeric, allowing discoveries [31], [32] to go further 
|33j . [33], [35]. Modern numerical relativity without cu- 
mulative theoretical insights, would not have developed 
successfully, of course. 

Among many possibilities of exchange between the two 
aforementioned approaches, some years ago were intro- 
duced seminumerical techniques, which may be regarded 
as a "compromise" between the analytical and numerical 
approaches. These techniques are based on a general al- 
gorithm for modeling self-gravitating spheres out of equi- 
librium and were initially developed for radiation (Bondi) 
non-comoving coordinates in [36] (see |37j for a review 
and further references). In that version the method has 
been applied to a variety of different physical scenarios 
(see for example [35] _ [1Z] and references therein). 

Later on the method was extended to Schwarzschild 
-like coordinates (non-comoving) in [48], [49]. In this 
format, the algorithm has also lead to a variety of ap- 
plications [5U]-[S3]. Technically, the original method, in 
radiation coordinates, is first order in the local radial ve- 
locity, whereas in Schwarzschild-like coordinates is sec- 
ond order. 

The proposed method (in either version), starting 
from any interior (analytical) static spherically symmet- 
ric ("seed") solution to Einstein equations, leads to a 
system of ordinary differential equations for quantities 
evaluated at the boundary surface of the fluid distribu- 
tion, whose solution (numerical), allows for modeling, dy- 
namic, self-gravitating spheres, whose static limit (when- 
ever it exists) is the original "seed" solution. 

The approach is based on the introduction of a set of 
conveniently defined "effective" variables (effective pres- 
sure and energy density) and an heuristic ansatzs on the 
later, whose rationale and justification become intelligi- 
ble within the context of the post-quasistatic approxi- 
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mation. In the quasistatic approximation, the effective 
variables coincide with the corresponding physical vari- 
ables (pressure and density) and therefore the method 
may be regarded as an iterative method with each con- 
secutive step corresponding to a stronger departure from 
equilibrium. 

It should be observed, that such seminumerical tech- 
niques require the possibility to approach the non- 
equilibrium state by means of succesive approximations, 
implying that there is life between quasi-equilibrium and 
non-equilibrium, at least in some cases. 

Motivated by the success of previous versions of the 
method in non-comoving coordinates, and by the fact 
that comoving coordinates are commonly used in the 
study of gravitational collapse, we endeavour in this work 
to look at some version of the above mentioned algorithm 
in comoving coordinates, which in turn would require the 
definition of the post-quasistatic approximation in co- 
moving coordinates. 

It concerns about the relationship between Eulerian 
and Lagrangian observers. From the historic point of 
view, is about Bondi's [M] and Misner-Sharp's [55] ap- 
proaches to deal with matter. The former leads to the 
Wilson codes, the latter to May and White ones [25J, in 
the context of modern numerical relativity |56j . 

In what follows we develop this plan: First, using co- 
moving coordinates we write down the field equations for 
the most general fluid. Second, we introduce appropri- 
ate definitions of mass and velocity, following Misner and 
Sharp, to recast the field equations. Third, we detail the 
junction conditions with the exterior spacetime, which 
is of Vaidya. Fourth, we consider the static, quasistatic 
and post-quasistatic regimes. Fifth, we propose a proce- 
dure for the modeling and illustrate the algorithm with 
a simple model based on the Schwarzschild interior so- 
lution. We consider the adiabatic and the nonadiabatic 
case. Finally we include some concluding remarks in the 
last section. 



where A, B and R are functions of t and r and are 
assumed positive. We number the coordinates x° = t, 
x 1 = r, x = 6 and x 3 = 4>. 

B. Energy— momentum tensor 

The matter energy-momentum T a p inside £ has the 
form 

T a p = (ji + P±)V a Vp + P±g a p + (P r -P±)XaXp 

+ q a Vp + V a qp + eljp - 2-qa af} , (2) 

where fi is the energy density, P r the radial pressure, 
P± the tangential pressure, q a the heat flux, e the en- 
ergy density of the null fluid describing dissipation in the 
free streaming approximation, rj the shear viscosity coef- 
ficient, V a the four velocity of the fluid, \ a a un h four 
vector along the radial direction and l a a radial null four 
vector. These quantities satisfy 

V a V a = -1, V a q a = 0, X a Xa = h 

x a v a = o, rv a = -l, n a = o. 

Observe that we have assumed the shear viscosity ten- 
sor 7r Q( 3 to satisfy the relation 

K a p = -2Wa/3, (3) 

where a a p is the shear tensor. However this last equa- 
tion is valid only within the context of the standard irre- 
versible thermodynamics (see [M], [55J for details). 

In a full causal picture of dissipative variables we 
should not assume Instead, we should use the trans- 
port equation derived from the corresponding theory (e.g. 
the Miillcr-Israel-Stewart theory [HH])- However for 
the sake of simplicity, in this work we shall restrict our- 
selves to the standard irreversible thermodynamics the- 
ory. 



II. COMOVING FRAMES TO DESCRIBE 
GRAVITATIONAL COLLAPSE 

A. Comoving coordinates 

We consider a spherically symmetric distribution of 
collapsing fluid, bounded by a spherical surface S. The 
fluid is assumed to be locally anisotropic (principal 
stresses unequal) and undergoing dissipation in the form 
of heat flow (to model dissipation in the diffusion approx- 
imation), null radiation (to model dissipation in the free 
streaming approximation) and shearing viscosity. Phys- 
ical arguments to consider such fluid distribution in the 
study of gravitational collapse may be found in [57]-[62] 
and references therein. 

Using comoving coordinates as in |63j . we write the 
line element in the form 

ds 2 = -A 2 dt 2 + B 2 dr 2 + R 2 (d9 2 + sin 2 9d(f> 2 ), (1) 



C. Kinematical variables 

The four-acceleration a a and the expansion of the 
fluid are given by 

a a = V a ., p V p , ® = V a , a , (4) 

and its shear o a p by 

Pafj = V(a;/3) + a (aVp) - -0/l a /3, (5) 

where h a p = g af} + V a Vp. 

We do not explicitly add bulk viscosity to the system 
because it can be absorbed into the radial and tangential 
pressures, P r and Pj_, of the collapsing fluid [70] , 

Since we assumed the metric ([l]) comoving then 

V a = A~ 1 8q, q a = qB- 1 ^, (6) 
l a = A~ l 5% +B- 1 5*, x a =B- 1 S?, (7) 
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where q is a function of t and r satisfying q a = qx° ■ 

From Q with ^ we have for the four-acceleration 
and its scalar a. 



A' 2 a M' 

ai= A> a =a aa ={jB 



(8) 



where 



°= l -A B - R 

A \ B R 



(12) 



where a a = ax a , and for the expansion 



Then, the shear tensor can be written as 



A\B R 



(9) 



where the prime stands for r differentiation and the dot 
stands for differentiation with respect to t. With Q we 
obtain for the shear ^ its non zero components 



Cq/3 = er ( XaX/3 - ^h a p 



(13) 



1 U 2 a 33 

o sm 6 



(10) 



D. Field equations 



and its scalar 



<7 o-q^ = -er 



Thus, the Einstein field equations for the interior 
(11) spacetime |l]) can be written as 



8TTflA 2 
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B R R \B 
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^B' R! 

'^B^R 



(14) 
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(15) 



8wP r B 2 = 
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(17) 



where 

M = M + e > 
4 

P r = P r - -7?(T + e, 

P± = P± + \w- 

Observe that if functions A(t,r), B(t,r) and R(t,r) are 
completely determined, the system above becomes an al- 
gebraic system of four equations for the six unknown 



functions [A, e, q, P r , P± and r/. In this general case ad- 
ditional equations are required (e.g. an equation of state 
and an equation describing energy production) in order 
to close the system. This is an expression of the well es- 
tablished fact that under a variety of circumstances a line 
element may satisfy the Einstein equations for different 
(physically meaningful) stress-energy tensors (see [71]- 
[52] and references therein). In the locally isotropic case 
dissipating in either the free streaming or the difussion 
limit the system is overdetermined. The same happens in 
the non-dissipative case, even if the fluid is anisotropic. 
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III. MASS AND VELOCITY 



Using ( 14 )-( 16 1 with ( 19 ) and ( 20 ) we obtain from ( 18 ) 



Following Misner and Sharp [55] , let us now introduce 
the mass function m(t,r) (see also [83]), defined by 



2 



p 23 



R 
2" 



R 



+ 1 



(18) 



It is useful to define the proper time derivative given 
by 



D T = 



Adt' 



and the derivative Dr, 



Dr 



l_d_ 
R/dr' 



(19) 



(20) 



where R defines the areal radius of a spherical surface 
inside £ (as measured from its area). 

Using ( 19 ) we can define the velocity U of the collaps- 



ing fluid as the variation of the areal radius with respect 
to proper time, i.e. 



U = D T R. 



Then ( 18 ) can be rewritten as 



2m\ 

~R) 



1/2 



(21) 



(22) 



D T m = -ATT P r U + qE R\ 



(23) 



and 



~U\ -2 



D R m = 4tt ( p, + q— ) R 



(24) 



Next, the three-acceleration DtU of an infalling par- 
ticle inside S can be obtained by using (16), (18) and 
( 22 ) , producing 



D T U = -£-4*P r R + E± i 



(25) 



t4/ 
~A 



4ttRB 
E 



D T U 



AttR 4ttR 3 



Pr 



(26) 



Now, from the Bianchi identities we obtain (see eq. 
(38) in 63J) in this case 



(p, + P r )D T U 



~ (P + Pr 



m 
R? 



4wP r R 



— E 



DoP r 



- E 



D T {e + q) + 2(e + q) I 



f2U 



(27) 



r 



The physical meaning of different terms in (27) has been 
discussed in detail in [5S]-[Sn]- Suffice to say in this point 
that the first term on the right hand side describes the 
gravitational force term. 



IV. THE EXTERIOR SPACETIME AND 
JUNCTION CONDITIONS 

Outside E we assume we have the Vaidya spacctime 
(i.e. we assume all outgoing radiation is massless), de- 
scribed by 



ds 2 



2M(v) 



1 dv 2 -2dpdv + p 2 (d6 2 + sin 2 6»<i^ 2 ), 

(28) 

where M(v) denotes the total mass, and v is the retarded 
time. 



The matching of the full nonadiabatic sphere (includ- 
ing viscosity) to the Vaidya spacetime, on the surface 
r = rs = constant, was discussed in [51] (for the discus- 
sion of the shear-free case see [SS] and 86J). However 
observe that we are now including a null fluid within the 
star configuration. 

Now, from the continuity of the first differential form 
it follows (see [84] for details), 



Adt = dv[l 



= dr, 



and 



2M(v) 



R k p(y), 



— \ = (l 2m 1 2 dp 
dr J \ p dv 



(29) 



(30) 



(31) 
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where r denotes the proper time measured on E. 



Whereas the continuity of the second differential form produces 



m(t, r) = M(v), 



and 



R 



BR' RA' 
~B~R ~ R~A 



e _B 
A 



'R 



A 



R\ R 
R I R 



r 



.4 
B 



A' R'\ R' 



A R R 



(32) 



(33) 



where = means that both sides of the equation are eval- 
uated on E (observe a misprint in eq.(40) in [84] and a 
slight difference in notation). 

Comparing (331 with (15) and (fl6| one obtains 



q = P r 



-r/a. 



(34) 



Thus the matching of ([lj and (28) on E implies (32) 
and (34), which reduces to equation (41) in [M] with the 
appropriate change in notation. Observe a misprint in 
equation (27) in [5U] (the a appearing there is the one 
defined in [M], which is —1/3 of the one used here and 
in |H]). 

Also, we have 



, e L 



(35) 



where Le denotes the total luminosity of the sphere as 
measured on its surface and is given by 



L = L n 



1 



2m 
P 



dp 
dv 



and where 



dM s 
dv 



dm dt . dv _ 1 
Ittdr^dr' 



(36) 



(37) 



is the total luminosity measured by an observer at rest 
at infinity. 

The boundary redshift z-s is given by 



dv s i , 
Jr =1 + Z > 



with 



dv E 
d7 = 



B 



R 
~A 



(38) 



(39) 



Therefore the time of formation of the black hole is given 
by 



B 



R 
A 



= E + U = i). 



(40) 



Also observe than from (31), (36) and (39) it follows 



(E + U f ' 



and from (21), (22), (31) and (39) 



dp e 
dv 



U(U + E). 



(41) 



(42) 



V. EVOLUTION REGIMES 

We shall next define three possible regimes of evolu- 
tion. 



A. Static regime 

In this case all time derivatives vanish, implying: 



U = 6 = a = 0. 



(43) 



Since B = B(r); A = A(r);R = R(r), reparametrizing 
r, we may write the line element in the form: 



ds 2 = -A 2 dt 2 



B 2 dr 2 



\d0 2 



sin 



; ). (44) 



Thus, the "Euler" equation ( |27| becomes the well 
known TOV equation of hydrostatic equilibrium for an 
anisotropic fluid 



PI 



-(P r P X ) = ( / i + f\ (m + 4nPy). (45) 
r r(r — 2m) 



The Einstein equations in this case read: 



8npA 2 = 



-2*- 
Br 



8TrP r B A 



i\ i 



^P,r 2 = Q' 



r \ 2 


'A" 


A 1 


B' 




B' 


b) 


A 


A 


B + 




B 



(46) 

(47) 
(48) 
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Also, for the mass function we have 
r ( 1 



B 2 = 1 - 



2m 



(49) 



(50) 



Thus in quasi-equilibrium we have to assume: 

0{U 2 ) = A 2 = B 2 = Ab = R = B^0 (53) 

and the radial dependence of the metric functions as well 
as that of physical variables is the same as in the satic 
case. The only difference with the latter case is the fact 
that variables depend upon time according to equation 
Jl5>. 



to = 4-7T / [ir dr, 
'a 



and for the metric function A, we have from (261 

(to + 47rr 3 P r ) 



In 



-dr. 



(51) 



(52) 



J Jr s r(r — 2m) 
Therefore, once the radial dependence of \i and P r are 



known, the metric functions arc determined from (50 
521. 



B. Quasistatic regime (QSR) 

As is well known, in this regime the system is assumed 
to evolve, but sufficiently slow, so that it can be consid- 
ered to be in equilibrium at each moment (Eq. ( 45 ) is 
satisfied). This means that the sphere changes slowly, 
on a time scale that is very long compared to the typ- 
ical time in which the sphere reacts to a slight pertur- 
bation of hydrostatic equilibrium, this typical time scale 
is called hydrostatic time scale [S7] - [HS] (sometimes this 
time scale is also referred to as dynamical time scale, 
e.g. [52]). Thus, in this regime the system is always very 
close to hydrostatic equilibrium and its evolution may be 
regarded as a sequence of static models linked by (15 1. 



This assumption is very sensible because the hydro- 
static time scale is very small for many phases of the life 
of the star [55]. It is of the order of 27 minutes for the 
Sun, 4.5 seconds for a white dwarf and 10 -4 seconds for 
a neutron star of one solar mass and 10 Km radius. It 
is well known that any of the stellar configurations men- 
tioned above, generally (but not always), change on a 
time scale that is very long compared to their respective 
hydrostatic time scales. Let us now translate this as- 
sumption in conditions to U and metric and kinematical 
functions. This implies that: 

• The areal velocity U as well as other kinematical 
variables are small, which in turn implies that dissi- 
pative variables and all first order time derivatives 
of metric functions are also small. 

• From the above and the fact that the system always 
satisfies the equation of hydrostatic equilibrium, it 
follows from p7| that second time derivatives of 
metric functions can be neglected. 



C. Post— quasistatic regime (PQSR) 

In the two regimes considered above the system is al- 
ways in (or very close to) hydrostatic equilibrium. Let 
us now move one step forward into non-equilibrium and 
let us assume then that ( 45 1 is not satisfied. 



Then the question arises: What is the closest situation 
to quasi-equilibrium, not satisfying eq. (451? For ob- 



vious reasons we shall call this regime, post-quasistatic 
regime. Three remarks are in order at this point: 

1. First of all it should be stressed that the main moti- 
vation to consider the PQSR is to have the possibil- 
ity to consider those aspects of the object directly 
related to the non-equilibrium situation, which for 
obvious reasons cannot be described within the 
QSR. 

2. It should be clear that we are also assuming the 
fact that we can approach the non-equilibrium by 
means of successive approximations. It goes with- 
out saying that not any self-gravitating fluid will 
satisfy this requirement. 

3. It also should be clear that unlike the two prece- 
dent regimes, there is not a unique definition for 
PQSR. In what follows we shall propose a defini- 
tion in analogy to the one given in [48 for non- 
comoving coordinates. 

4. Once the system is out of equilibrium, its eventual 
return to the static or quasistatic regime is not as- 
sured and will depend on the initial data and the 
very nature of the system. 

Now, since in both, the static and quasistatic regimes, 
the radial dependence of metric variables is the same, we 
shall keep that radial dependence as much as possible, 
but of course the time dependence of those variables is 
such that now ( 53 ) is not satisfied. 



Then from the above we write 



R = rn(t), 



(54) 



where k is an arbitrary function of t, to be determined 
later. 

Taking into account (22) and (54), we rewrite the met- 
ric as follows 



ds 2 



A 2 dt 2 + K 2 [E~ 2 dr 2 + r 2 (d6 2 + sin^d^)]. (55) 
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Next, defining the effective mass as 

1 



we obtain 



l eff 



E 2 



RU 2 , 



1 



2m, 



If 



R 



Then, equations ( 24 1 and ( 26 ) can be written as 
1 



-m 



'eff = 4nR Veff, 



with 



-(In A)' 

K 



Veff 



Peff = Pr 



4irR 2 P ( 



eff 



eff/R 



R-2m 



eff 



qU _ UD R U 
E 



U 2 



AttR 8ttR 2 ' 

D T U U 2 
AttR + 8ttR 2 ' 



(56) 
(57) 

(58) 
(59) 

(60) 
(61) 



where we have followed the terminology used in [36] , [37] , 
[35] and call /tt e // an d Peff the "effective density" and 
the "effective pressure", respectively. The meaning of 
these variables will become clear in the discussion below. 

It is worth observing that in comoving coordinates the 
PQSR is "exactly" second order in U (or q), which is 
not the case in non-comoving coordinates (radiation or 
otherwise) . 



Next, from ( 58 )— ( 61 ) , with (54) we may write 



1 f r , 

-Z3 m eff = / 47!T Veffdr, 
K Jo 



(62) 



^R 3 Peff + rn e ff) 

R(R-2m eff ) 



dr. (63) 



n(t)r 



From the above, it follows at once that if R 
and [J-eff shares the same radial dependence as /i in the 
static case, then obviously the radial dependence of m e // 
will be the same as in the static case. The inverse is true 
of course, if the radial dependence of rn e ff is the same 
as in the static case, then fJ. e ff shares the same radial 
dependence as \i static. 

On the other hand, if besides the assumption above, 
we assume that P e ff shares the same radial dependence 
as P r static, then it follows from (631 that A shares the 



same radial dependence as in the static case. 

All these considerations provide the rationale for our 
algorithm. Indeed, starting with a "seed" static solution 
with a given n(r) and P r (r), it follows from (62 1 and (63 1, 
that the PQSR can be implemented by assuming that the 
radial dependence of the effective variables is the same 
as that of /x(r) and P r {r) (thogether with R = rn(t)). 



It is worth stressing an important difference between 
the operational definition of QSR and PQSR in comov- 
ing and noncomoving frames. In the latter case there is 
one physical variable more (the velocity) and one metric 
function less (i?) than in the former. 

We shall next outline the approach that we propose. 

As mentioned before, such an approach was already 
proposed and developed in a noncomoving frame of refer- 
ence (see [36] , [37] and [48] and references therein) . Here 
we want to provide a formulation for comoving frame. 

The proposed method, starting from any interior (ana- 
lytical) static spherically symmetric ( "seed" ) solution to 
Einstein equations, leads to a system of ordinary differ- 
ential equations for quantities evaluated at the bound- 
ary surface of the fluid distribution, whose solution (nu- 
merical), allows for modeling, dynamic, self-gravitating 
spheres, whose static limit is the original "seed" solution. 

The approach is based on the post-quasistatic assump- 
tion, which as mentioned before is equivalent to the 
assumption that "effective" variables (effective pressure 
and energy density) share the same radial dependence 
as the radial pressure and energy density of the static 
"seed" solution, respectively. An ansatzs justified by the 
fact that in the quasistatic approximation, the effective 
variables coincide with the corresponding physical vari- 
ables (pressure and density). Therefore the method may 
be regarded as an iterative method with each consecutive 
step corresponding to a stronger departure from equilib- 




FIG. 1: Evolution of the radius i?s for model A and initial 
conditions Jfe(0) = 10; m s (0) = 1; (7e(0) = 0; ^(0) = 
-10" 2 . 
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FIG. 2: Evolution of the energy density fi = m E (0)/i (mul- 
tiplied by 10 2 ) for model A and initial conditions Rs(0) = 

10; m E (0) = 1; £fc(0) = 0; ^(0) = -1(T 2 . 



FIG. 4: Evolution of the radius _R E for model C and initial 

conditions ife(0) = 10; m E (0) = 1; t/ E (0) = 0; ^(0) = 

dr 

-10" 2 . 




-0.3 - 

-0.4 - 

-0.5 - 

-0.6 - 

-0.7 - 

-0.8 - 

-0.9 - 



r/r 2 =0 


2 


r/r 2 =0 


4 


r/r 2 =0 


6 


r/r 2 =0 


8 


r/r y =l 










10 



15 20 

X 



25 30 



35 




FIG. 3: Evolution of the velocity U for model A and initial 

conditions ife(0) = 10; m E (0) = 1; Us(0) = 0; ^(0) = 

dr 

-10" 2 . 



FIG. 5: Evolution of the total mass m E for model C and initial 
conditions Jfe(0) = 10; m E (0) = 1; t/ E (0) = 0; ^(0) = 
-10~ 2 . 



VI. PROTOCOL 

On the basis of all comments above, we shall now 
present the following algorithm: 

1. Take an interior ( '"seed" ) solution to Einstein equa- 
tions, representing a fluid distribution of matter in 
equilibrium, with a given 



Assume that the r dependence of effective variables 
is the same as that of P r and \i respectively and 
R = rn(t). This assures that the remaining metric 
functions share the same radial dependence as that 
of the "seed" solution. 



Using equations ( 62 1 and ( 63 ) , with the r depen- 
dence of P e ff and Heff, one gets m e ft and A up 
to some functions of t. 



H = fi(r); P r = P r (r). 



4. For these functions of t one has three ordinary dif- 
ferential equations (hereafter referred to as Surface 
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FIG. 6: Evolution of U-z[dUv/df]_ (multiplied by 10 4 ) for 
model C and initial conditions 7?e(0) = 10; ms(0) = 

1; £/ E (0)=0; ^(0) = -1CT 2 . 



equations), namely: 



(a) Equation (21) evaluated on r = ry. 

(b) Equation ( 23 ) evaluated on r = rs . 

(c) Equation (27) evaluated on r = ?*£. 



5. Depending on the kind of matter under consider- 
ation, the system of surface equations described 
above may be closed with the additional informa- 
tion provided by the transport equation and/or 
the equation of state for the anisotropic pressure 
and/or eventual additional information about some 
of the physical variables evaluated on the boundary 
surface (e.g. the luminosity). 

6. Once the system of Surface equation is closed, it 
may be integrated for any particular initial data. 

7. From the result of integration we obtain R and then 
using (22), (62) and (63), m e ff, B and A, are com- 
pletely aetermined. 

8. With the input from the point 7 above, and us- 
ing field equations, together with the equations of 
state and/or transport equation, all physical vari- 
ables may be found for any piece of matter distri- 
bution. If the system is "very far" from equilibrium 
then it could be necessary to go through the pro- 
cess once again by replacing the seed solution by 
the solution obtained from the point 7 above. This 
could be done as many times as it is required to ob- 
tain a satisfactory description of the system under 
consideration. 

Let us now elaborate on the surface equations, which 
are the corner stones of the proposed algorithm. 



As mentioned above, the first surface equation is equa- 

(64) 



tion (21 ) evaluated on r = r%, i.e. 

Ry = U s A? 



or using ( 29 ) 



dr 



U, 



(65) 



It would be convenient to introduce the dimensionless 
variables 

t 



t = 



m s (0)' 

T 

m s (0)' 
= m E (0) : 



and 



Ry, = 



m E (0)' 



(66) 
(67) 
(68) 

(69) 



where ms(0) denotes the total initial mass. 
Then, the first surface equation reads 



or 



dRy. 
di 



dRy 
df 



UyAs 



Uy. 



(70) 



(71) 



Next, evaluating (23) on £ and using (34), we get 



D T m = -4tt [(q + e)(U + E)] R 2 



or using ( 30 ) and ( 35 1 



D T m = -(U + E)L. 



(72) 



(73) 



In terms of dimensionless variables this last equation 
reads 



dm^ 
~d? 

drily: 
di 



= -(U + E)vLv, 
-- -A S (U + E) s Ls. 



(74) 
(75) 



Equation (74) (or ( [75] )) is the second surface equation. 
Instead of working with Ly (the luminosity of the object 
as measured on its surface), we may replace it by the 
luminosity measured by an observer at infinity , using 
JTT:, 



The third surface equation may be obtained from the 
"Euler equation" (27) evaluated at the boundary surface. 



The general form of this equation is quite long, and we 
shall only write it explicitly for the simplified models con- 
sidered below. 



10 



2.75 r 

2.70 - 

2.65 - 

2.60 - 

1^-2.55 - 

2.50 - 

2.45 - 
2.40 
2.35 



r/r 2 =0 


2 


r/r 2 =0 


4 


r/r 2 =0 


6 


r/r 2 =0 


8 


r/r y =l 


- 





FIG. 7: Evolution of the energy density ft = m E (0)/i (mul- 
tiplied by 10 4 ) for model C and initial conditions -R E (0) = 

-2 



10; rn E (0) = 1; C7 S (0) = 0; 
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FIG. 9: Evolution of the grade of anisotropy A = P± — P r = 
■m%(0)(P± ~ -Pr) (multiplied by 10 6 ) for model C and initial 

conditions ife(0) = 10; m E (0) = 1; U s (0) = 0; ^(0) = 

-10" 2 . 
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FIG. 8: Evolution of the radial pressure P r = m E (0)P r (mul- 
tiplied by 10 ) for model C and initial conditions -Rs(O) = 

10; r?te(0) = 1; (7e(0) = 0; ^(0) = -10~ 2 . 



FIG. 10: Evolution of the energy flux e = m s (0)e (multi- 
plied by 10 6 ) for model C and initial conditions -R E (0) = 

10; m E (0) = 1; E7 E (0) = 0; ^(0) = -10"" 2 . 

UT 



VII. MODELING: SOME EXAMPLES 

We shall now illustrate the method outlined above with 
some examples inspired in the well known Schwarzschild 
interior solution. It should be clear that our goal here is 
not to solve any specific astrophysical problem, but just 
to exhibit the potential of the method. Also it should 
be mentioned that models inspired in the Schwarzschild 
interior solution (within the context of the PQSR) were 
presented in [35], [35] and [25] , using the version of the 



method in non-comoving coordinates. However, the ini- 
tial data as well as the type of anisotropy used here differ 
from the one assumed for those previous models, and ac- 
cordingly the evolution pattern is quite different. This 
latter point strees further the fact that in spite of the 
similarity of the basic assumption underlying the very 
definition of the PQSR, the framework of both versions 
(in comoving and non-comoving frames) are are very dif- 
ferent. 

For the sake of simplicity we shall consider nonviscous 
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FIG. 11: Evolution of the velocity U for model C and initial 

conditions ife(0) = 10; m E (0) = 1; !7 S (0) = 0; -t#(0) = 

ar 

-10" 2 . 



fluid, dissipating only in the streaming out limit (i.e. q = 

n = o). 

Thus, the effective density for the Schwarzshild type 
model is assumed to depend only on the time-like coor- 
dinate, i.e. 

Me// = /(*), (76) 
whereas the expression for the effective pressure is 

p eff + 



(J-eff + Peff 



= £A(*), 



with 



1/2 



We may also write ( 77 1 as 



P eff = faff 



3(X - &) ' 



(77) 



(78) 



(79) 



where A = 4>/x is an arbitrary function of time such that 
the pressure satisfy the boundary condition, which in this 

case is P r = 0. 

It is easy to check that 



ar 



X = E^{ i? s L s + 3m s - RvU£ 



dU- 



df 



s m 



(80) 



(81) 



Also, using (62) and (63) we obtain 

m eff = y^ 3 /' 



A = Ay 



where 



and 



Then (27) is written as 



/2ms _ 
I -Re 



U4 



(82) 



(83) 



(84) 



(85) 



(fi + P r )D T U + (fi + P r ) UwRP r 



m \ 
R?) 



E D R P r 



+E 



D T e + 2e 



R P - 
2U 
~R~ 



(86) 



Next, we shall evaluate (86) at the boundary surface, 



for doing that we shall need an equation of state for the 
stresses (at the boundary). For the sake of simplicity 
we consider [Pj_]e = 0, then after lenghty manipulations 
using MAPLE, we are lead straightforwardly to the third 
equation at the surface 



JT d 2 Uv . fdUj; 
Us ,_o = ■A 



where 



df 2 



A = 



df 



(87) 



J_ / 2m s 
El \ Rz 



-2UZ-1 



B 



(-3C/|i?| - 3I/p| + llt/|m s i?| 



+ 2m|i? s - m s i?| + [D T m]TjJnRl) > 



(89) 



1 



[L s (2i?| + 2(7|i?| 



8m| - 8C/p s m E 



R% E s 

+4f/p| - 8i? s rn s ) + [£> T m] E (8™ s i? s [/ s 



-RlU. 



E u| - i?|(7 s - 3[D r m] E J^j) 



(90) 



The system of surface equations (71), (74) and (87) 



has been integrated for the following initial conditions 
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i? E (0) = 10; m s (0) = 1; C/ s (0) = 0; ^(0) = -HT 2 

(XT 

We shall consider both, the adiabatic and the dissi- 
pative case (in the free streaming approximation), in the 
latter case we shall use a Gaussian as a luminosity profile 



Ly 



L e~ 



We integrate numerically the system of equations at the 
surface using a standard fourth order Runge-Kutta. We 
built up a nonadiabatic model preparing as test beds two 
previous models do not deprived of physical interest. 

For simplicity we have cut here the iterative process 
after the first step of the chain (i.e. we have assumed 
the seed solution as given by (76 1 and (79)). We invite 
any interested reader to proceed to the next step of the 
iterative method by replacing the effective variables by 
the energy density and pressure obtained here in the first 
step and go through the protocol once again. 



A. Model A: Fully adiabatic 

In this model we use DtB = BDrU everywhere which 
is equivalent to enforce e = 0. Thus, setting Lq = wc 
obtain results displayed in figures 1-3. The isotropic fluid 
sphere behaves as dust, that is P r = P± = everywhere, 
and collapse proceeds catastrophically as expected. The 
monitoring of the evolution is stopped when trespassing 
the horizon. 



B. Model B: Without luminosity 

We can set only Lq = to see how the distribution 
evolves. With no surprise the results are the same as in 
model A. This is a consequence of the well-posed initial- 
boundary problem. 



C. Model C: Nonadiabatic 

When some fraction of the total mass (« 1%) is car- 
ried away by the Gaussian pulse, picked at fo = 5, with 
£ = 1, the dust ball initially goes to collapse, becoming 



anisotropic in the process. These results are shown in 
figures 4-11. 



VIII. CONCLUDING REMARKS 

A seminumerical method to describe gravitational col- 
lapse in comoving coordinates has been proposed, in anal- 
ogy with the already existing algorithm in non-comoving 
coordinates. 

For doing that, we have revisited and redefined the 
basic concepts of the post-quasistatic approximation 
(PQSA) in order to adapt them to comoving coordinates. 
The essential features of the seminumeric method keep 
going. But in comoving coordinates we were enforced 
to reasonable transfer the PQSA to a geometrical vari- 
able (R). Up to now the effective variables let us to 
make heuristically the job. In the present version an ad- 
ditional geometrical point of view led straightforwardly 
to the new effective variables. We endeavor supposing 
there is life in between quasistatic and post-quasistatic 
regimes. So we did. Here we reported one version of the 
PQSA in comoving coordinates. 

Wc have integrated the surface equations for some sim- 
ple models inspired in the interior Schwarzschild solution. 
Our intention presenting such models was not to describe 
any physically relevant astrophysical scenario, but just to 
illustrate the method. Due to the simplifications imposed 
on the models, the field equations are overdetermined 
producing specific constraints. Accordingly, a fine tun- 
ing specification of initial conditions is necessary. Even 
more, in our models we never recover the quasistatic o 
the static regimes. Our models are intrinsically unsta- 
ble and physically acceptable. We have to mention that 
we observe certain tendency to stabilize the system if the 
configuration initially was less relativistic (less compact). 

Emission of energy seems to play a crucial role in the 
process of gravitational collapse. Eventually dissipation 
and anisotropy (unequal stresses) may change the evolu- 
tion fate, avoiding the complete collapse to a black hole 
in the same hydrodynamical time scale. Some additional 
work is required on this last issue, specially considering 
a transport equation to deal with heat flow (and/or vis- 
cosity) dissipation in the context of extended thermody- 
namics within the PQSA (see [27], [HO] for a treatment 
of this problem in non-comoving coordinates). 
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